Housing market price estimators like Zillow’s Zestimate can give a general idea of a house’s value, but more often than not, professional appraisers and real estate agents do no recommend relying on these tools as they tend to have high error rates. In addition to inconveniencing a home buyer or seller, the high error rates of automated valuation models (AVM; e.g. Zillow’s Zestimate) may misrepresent low-income and majority-minority, especially Black-majority neighborhoods where home price value errors are higher compared to majority-white neighborhoods1. As well, a home’s value is often used to determine a homeowner’s mortgage loan approval and to inform policymakers on the general state of the housing market 2. Inaccurate housing estimates may lead to misinformed housing and economic policies disproportionately affecting non-white neighborhoods, reifying racial and class disparities in housing.
There are many reasons why these price estimators may be inaccurate (e.g. data availability, inability to capture temporal housing market patterns, etc.), and among them lack of local intelligence stands out. Each urban area has its own culture and unique characteristics that need to be factored into AVMs in order to more accurately predict housing prices. As residents and spatial analysts of the city of Philadelphia, our team was able to apply our knowledge of the city and its neighborhoods to the creation of an improved prediction model.
Even with a deeper understanding of Philadelphia, creating an improved home price estimator model for the city is still a challenging task. All kinds of data are available including health and human services, transportation, public safety, parks and recreation, and food in addition to the traditional housing market and attribute data used in AVMs. Though all of these data are interrelated and inform each other and home price, limiting the amount of predictors to only keep those that are most informative can be a daunting task. We often create new features by redefining (sometimes by combining) variables that improve our prediction of home price. Luckily we implement technology that makes this decision-making process easier.
Another possible challenge when it comes to constructing a hone price estimator model is data availability, or rather unavailability. Some neighborhoods, especially low-income neighborhoods and regions where the population does not trust in government, data that we use in our AVMs may be missing for various reasons. Some of the variables that are most informative to a model may be resources that are unavailable in low-income neighborhoods. On the other hand, even if the resources are available and utilized in the community, a region with high distrust in government may not report information in the government Census and other data surveys, resulting in missing data. Ultimately, if a model relies on variables that are missing in neighborhoods with certain characteristics, it will predict home price with more error in these regions.
Keeping these potential pitfalls in mind, our overall modeling strategy was as follows. We first gathered open-source data about Philadelphia to examine and create variables of Philadelphia houses and neighborhoods to assess whethere it made sense to include them in our model. We then input a selection of predictors into an ordinary least squares (OLS) regression, specifically an OLS step-wise regression, to narrow down and determine the combination of the most informative predictors. After assessing our model fit by looking at mean absolute error (MAE) and mean absolute percentage error (MAPE), we revised our selection of predictors and repeated the process until we were happy with our model.
Our newly improved model uses 21 independent variables to predict housing prices in Philadelphia and is able to do so with a mean MAE of $68,687.82 and mean MAPE of 36.74% (with 100-fold cross-validation; see Table 7). Notable predictors in our OLS regression are number of bathrooms, percent of population that is white (in the census tract that the house is part of), and the rating of a house’s interior condition (see table 5). We observed clustering in home price predictions with a Moran’s I of 0.615 (p < 0.001) meaning our model errors exhibit spatial autocorrelation (see figure 10). We found that our model predicts with higher errors for majority non-White neighborhoods and low-income neighborhoods (see tables 8 and 9).
We gathered data about the shape and outline of the city of Philadelphia and its neighborhoods from OpenData arcGIS and Azavea’s Github repository3, respectively. We also pulled information about commercial centrs from OpenData arcGIS. Housing market and attribute data (including house features, sale price, ownership, etc.) was provided by Zillow (aka studentData), and additional census tract data (including total population, median income, etc.) was gathered from 5-year ACS Census data for 2020. We would like to note here that we previously had gathered data from the American forests website and Philly Carto websites about tree canopy cover and shootings, respectively, but our exploratory data anlalysis found that these data did not help us improve our model.
# data %>%# group_by(mapname) %>%# summarize(avg_sale_price = mean(sale_price)) %>%# arrange(desc(avg_sale_price)) %>%# head()# m <- cor(numeric_only %>% na.omit())# corrplot::corrplot(m, na.action = na.exclude)# numeric_only <- data %>% st_drop_geometry() %>% select(where(is.numeric))# ols_step_best_subset(model)# model <- lm(sale_price ~ ., data = numeric_only)# aic_preds <- ols_step_both_aic(model)$predictorskeep_vars <-c("price_lag","garage_spaces","total_livable_area","medHHInc","number_of_rooms","number_of_bathrooms", "depth","off_street_open","fireplaces","interior_condition","number_stories","pctWhite","number_of_bedrooms","exterior_condition","dist_to_commerce","year_built","totPop","total_area","sale_price", "hood_price_class", "building_type_price_class", "quality_grade","musaID","toPredict") # keeping toPredict to split data after creating dummy columnsdata_to_dummy <- data %>%select(all_of(keep_vars)) # for making summary statistics, keep categorical data as is and add imputed values# count sum of missing values# na_count <- c()# for (i in 1:ncol(data_to_dummy)) {# row1 <- data_to_dummy[,i]# na_count <- c(na_count,sum(is.na(row1)))# }# # na_count <- data.frame(name = names(data_to_dummy), na_count) %>% # filter(na_count > 0)# columns with missing data: garage_spaces, medHHInc, number_of_rooms, number_of_bathrooms, depth, off_street_open, fireplaces, interior_condition, number_stories, number_of_bedrooms, exterior_condition, total_areadummied_data <- data_to_dummy %>%dummy_cols(select_columns =c("hood_price_class","building_type_price_class","quality_grade")) %>%clean_names() %>%select(-c(hood_price_class, building_type_price_class, quality_grade))# impute missing data# to do this, have to remove to_predict (no characters are allowed + this has no NAs)# convert binary + integer (number of rooms) data into factors before imputing# there are some rows with 0 or 1 for total_area and total_livable_area. since this is impossible, we're replacing them with NA and imputing# dummied_data_toimpute <- dummied_data %>% # st_drop_geometry() %>% # dplyr::select(-c(geometry,to_predict)) %>% # mutate(total_area = ifelse(total_area < 2, NA, total_area),# total_livable_area = ifelse(total_livable_area < 2, NA, total_livable_area),# garage_spaces = as.factor(garage_spaces),# number_of_rooms = as.factor(number_of_rooms),# number_of_bathrooms = as.factor(number_of_bathrooms),# number_of_bedrooms = as.factor(number_of_bedrooms),# fireplaces = as.factor(fireplaces),# interior_condition = as.factor(interior_condition),# exterior_condition = as.factor(exterior_condition),# number_stories = as.factor(number_stories))# method without turning num vars to factors, round to nearest integer where necessary after imputation# number of bedrooms, bathrooms, and rooms == 0 need to be replaced with avg values# year_built == 0 needs to be replaced with average yearavg_rooms <- dummied_data %>%filter(number_of_rooms >0, !is.na(number_of_rooms)) %>% dplyr::select(number_of_rooms) %>%summarise(avg_rooms =mean(number_of_rooms)) %>%pull(avg_rooms) %>%round()avg_beds <- dummied_data %>%filter(number_of_bedrooms >0, !is.na(number_of_bedrooms)) %>% dplyr::select(number_of_bedrooms) %>%summarise(avg_bedrooms =mean(number_of_bedrooms)) %>%pull(avg_bedrooms) %>%round()avg_bath <- dummied_data %>%filter(number_of_bathrooms >0, !is.na(number_of_bathrooms)) %>% dplyr::select(number_of_bathrooms) %>%summarise(avg_bathrooms =mean(number_of_bathrooms)) %>%pull(avg_bathrooms) %>%round()avg_year <- dummied_data %>%filter(year_built >0, !is.na(year_built)) %>% dplyr::select(year_built) %>%summarise(avg_year =mean(year_built)) %>%pull(avg_year) %>%round()dummied_data_toimpute2 <- dummied_data %>%st_drop_geometry() %>% dplyr::select(-c(geometry,to_predict)) %>%mutate(number_of_rooms =ifelse(is.na(number_of_rooms) | number_of_rooms ==0,avg_rooms,number_of_rooms),number_of_bedrooms =ifelse(is.na(number_of_bedrooms) | number_of_bedrooms ==0,avg_beds,number_of_bedrooms),number_of_bathrooms =ifelse(is.na(number_of_bathrooms) | number_of_bathrooms ==0,avg_bath,number_of_bathrooms),year_built =ifelse(is.na(year_built) | year_built ==0,avg_year,year_built),total_area =ifelse(total_area <2, NA, total_area),total_livable_area =ifelse(total_livable_area <2, NA, total_livable_area))# imp_data <- missForest(dummied_data_toimpute,ntree = 20)$ximp %>%# cbind(dummied_data %>% dplyr::select(to_predict))# when rendering, suppressing this chunk for now for sake of time# imp_data2 <- missForest(dummied_data_toimpute2,ntree = 20)$ximp %>%# cbind(dummied_data %>% dplyr::select(to_predict))# saving imp_data so I don't have to run it again# write_csv(imp_data, "imputed_data_wfactors_231010.csv")# write_csv(imp_data2, "imputed_data_231011.csv")# imp_data <- read_csv("imputed_data_wfactors_231010.csv")imp_data2 <-read_csv("imputed_data_231011.csv")# for imp_data2, round to nearest whole number for garage_spaces, number_of_rooms, number_of_bathrooms, fireplaces, interior_condition,number_stories,number_of_bedrooms,exterior_conditionimp_data <- imp_data2 %>%mutate(garage_spaces =round(garage_spaces),number_of_rooms =round(number_of_rooms),number_of_bathrooms =round(number_of_bathrooms),fireplaces =round(fireplaces),interior_condition =round(interior_condition),number_stories =round(number_stories),number_of_bedrooms =round(number_of_bedrooms),exterior_condition =round(exterior_condition))reg_data <- imp_data %>%filter(to_predict =="MODELLING") %>%# only keep modelling dataselect(-c(to_predict))# geom and musaID to add back laterreg_data_geom <- data_to_dummy %>%filter(toPredict =="MODELLING") %>% dplyr::select(musaID,geometry)
Exploratory Analysis
We first limited studentData to rows with realistic values for sale price (sale_price > 1). After combining all data described above (matching by geographical location), we explored the relationship between variables and sale price. Scatterplots (like those of Figure 1) plot sale price as a function of interesting predictors. Plots B and D show that total livable area (measured in square ft) and the spatial lag of home price (using 15 nearest neighbors) have strong positive relationships with sale price. As expected, as median household income (measured at the census tract level) increases, so does the sale price of a home, though this positive relationship is not as strong as that with the previously mentioned features. Another interesting relationship we see in figure 1 is that of distance to nearest commercial center and sale price. Though one might expect a strong positive correlation between these two variables, the data shows a rather weak relationship between the two.
Show the code
# below is just a template of the codest_drop_geometry(reg_data) %>% dplyr::select(sale_price,dist_to_commerce,price_lag,total_livable_area,med_hh_inc) %>%rename(`A. Distance to Commercial Center (ft)`= dist_to_commerce,`B. Total Livable Area (sq ft)`= total_livable_area,`C. Median Household Income ($)`= med_hh_inc,`D. Spatial Lag of Home Sale Price`= price_lag) %>%gather(Variable, Value, -sale_price) %>%ggplot(aes(Value, sale_price)) + ggplot2::geom_point(size = .5, color ="#6e2e10", alpha =0.4) +geom_smooth(method ="lm", se=F, colour ="#027a59") +facet_wrap(~Variable, ncol =2, scales ="free") +labs(title ="Price as a function of 4 variables of interest",y ="Home Sale Price", x ="",caption ="Figure 1") +theme_minimal()
Feature Engineering
We observed some impossible or missing values in our data, namely values of zero or one in total_area and total_livable_area, missing values for spatial lag in price, zero for number of rooms, bedrooms and bathrooms, and zero for year_built. Though we knew why some of these variables were missing (e.g. spatial lag of price was missing for those that had a house in the challenge set as one of the 15 nearest neighbors since we set sale price for those houses as missing for this calculation), we had to get creative for the other cases. The different solutions we implemented were replacing impossible values as NA and imputing more reasonable values (in the case of total_area and total_livable_area) and replacing zeros with the mean value of that variable (for year_built, number of rooms, number of bathrooms, and number of bedrooms).
For categorical variables such as building type price class, quality grade, and neighborhood price class, we created dummy variables. Dummy variables help us break down each category of a categorical variable into its own feature, usually made up of ones and zeros indicating whether that row has that feature or not, respectively.
We have provided summary and descriptive statistics of the final predictors we decided to include in our model (tables 1-4).
Show the code
# reorder numeric datasum_stat_data <- data_to_dummy %>%st_drop_geometry() %>% dplyr::select(-c(garage_spaces, medHHInc, number_of_rooms, number_of_bathrooms, depth, off_street_open, fireplaces, musaID, interior_condition, number_stories, number_of_bedrooms, exterior_condition, total_area, sale_price, price_lag, total_livable_area, year_built)) %>%cbind(imp_data %>% dplyr::select(garage_spaces, med_hh_inc, number_of_rooms, number_of_bathrooms, depth,price_lag, off_street_open, fireplaces,interior_condition, number_stories,number_of_bedrooms, exterior_condition, total_area,total_livable_area, year_built)) %>% dplyr::select( # select in the desired order total_area, total_livable_area, number_of_rooms, number_of_bathrooms, number_of_bedrooms, number_stories, garage_spaces, fireplaces, depth, interior_condition, exterior_condition, year_built, quality_grade, building_type_price_class, off_street_open, # internal characteristics dist_to_commerce, # amenities/public services price_lag, med_hh_inc, pctWhite, totPop, hood_price_class # spatial variables )# summary statistics for numeric variablessum_stat_table <- sum_stat_data %>% dplyr::select(-c(building_type_price_class,quality_grade,hood_price_class)) %>% psych::describe() %>% dplyr::select(mean:median,min,max) %>%mutate(Mean =as.character(round(mean)), SD =as.character(round(sd)), Median =as.character(round(median)), Min =as.character(min), Max =as.character(round(max))) %>% dplyr::select(-(mean:max)) %>%filter(Max >1) %>%as.data.frame()# keep pctWhite rounded to nearest thirdsum_stat_pctwhite <- sum_stat_data %>% dplyr::select(pctWhite) %>% psych::describe() %>% dplyr::select(mean:median,min,max) %>%mutate(Mean =as.character(round(mean,3)), SD =as.character(round(sd,3)), Median =as.character(round(median,3)), Min =as.character(min), Max =as.character(round(max,3))) %>% dplyr::select(-(mean:max)) %>%as.data.frame()sum_stat_table <-bind_rows(sum_stat_table[1:16,], sum_stat_pctwhite, sum_stat_table[17,])sum_stat_table$Description <-c("Total area in square feet","Total liveable area in square feet","", "", "", "","Number of garage spaces","Number of fireplaces","Depth, as measured from the principal street back to the rear property line or secondary street","Overall condition of the interior","Overall condition of the exterior","", "","Distance to nearest commercial center","Spatial lag of sale price, 15 nearest neighbors","Median household income, census tract level data","Percent of population that is white, census tract level data","Total population, census tract level data")sum_stat_table %>%kbl(caption ="Table 1: Summary Statistics of numeric variables", align =rep("c", 8)) %>%kable_classic(full_width = F, html_font ="Cambria")
Table 1: Summary Statistics of numeric variables
Mean
SD
Median
Min
Max
Description
total_area
1920
3795
1241
150
226512
Total area in square feet
total_livable_area
1336
547
1210
64
15246
Total liveable area in square feet
number_of_rooms
4
1
4
1
32
number_of_bathrooms
1
1
1
1
8
number_of_bedrooms
3
1
3
1
31
number_stories
2
1
2
1
5
garage_spaces
0
1
0
0
72
Number of garage spaces
fireplaces
0
0
0
0
5
Number of fireplaces
depth
76
36
73
0
1115
Depth, as measured from the principal street back to the rear property line or secondary street
interior_condition
4
1
4
0
7
Overall condition of the interior
exterior_condition
4
1
4
1
7
Overall condition of the exterior
year_built
1938
30
1925
1750
2024
off_street_open
2005
1705
1535
14
12685
dist_to_commerce
670
638
546
0
7384
Distance to nearest commercial center
price_lag
268958
174639
226359
35440
1913110
Spatial lag of sale price, 15 nearest neighbors
med_hh_inc
54550
25387
50227
9276
155000
Median household income, census tract level data
pctWhite
0.438
0.315
0.428
0
1
Percent of population that is white, census tract level data
totPop
4691
1690
4473
0
9146
Total population, census tract level data
Show the code
# descriptive statistics for categorical datasum_stat_data %>% dplyr::select(building_type_price_class) %>%mutate(building_type_price_class =factor(building_type_price_class, levels =c("least expensive","less expensive","expensive","more expensive","most expensive"))) %>%table() %>%as.data.frame() %>%mutate(Description =c("Building types with average sale price of under 200,000 USD","Building types with average sale price between 200,000.01 and 400,000 USD","Building types with average sale price between 400,000.01 and 600,000 USD","Building types with average sale price between 600,000.01 and 800,000 USD","Building types with average sale price of over 800,000.01 USD")) %>%kbl(caption ="Table 2: Descriptive statistics of sales price classes of building type", align =rep("c", 8), col.names =c("Building Type", "Row count", "Description")) %>%kable_classic(full_width = F, html_font ="Cambria")
Table 2: Descriptive statistics of sales price classes of building type
Building Type
Row count
Description
least expensive
6679
Building types with average sale price of under 200,000 USD
less expensive
15484
Building types with average sale price between 200,000.01 and 400,000 USD
expensive
1156
Building types with average sale price between 400,000.01 and 600,000 USD
more expensive
552
Building types with average sale price between 600,000.01 and 800,000 USD
most expensive
10
Building types with average sale price of over 800,000.01 USD
Show the code
sum_stat_data %>% dplyr::select(hood_price_class) %>%mutate(hood_price_class =factor(hood_price_class, levels =c("least expensive","less expensive","expensive","more expensive","most expensive"))) %>%table() %>%as.data.frame() %>%mutate(Description =c("Houses in a neighborhood with average sale price of under 100,000 USD","Houses in a neighborhood with average sale price between 100,000.01 and 250,000 USD","Houses in a neighborhood with average sale price between 240,000.01 and 500,000 USD","Houses in a neighborhood with average sale price between 500,000.01 and 750,000 USD","Houses in a neighborhood with average sale price of over 750,000.01 USD")) %>%kbl(caption ="Table 3: Descriptive statistics of sales price classes of nieghborhoods", align =rep("c", 8), col.names =c("Neighborhood", "Row count", "Description")) %>%kable_classic(full_width = F, html_font ="Cambria")
Table 3: Descriptive statistics of sales price classes of nieghborhoods
Neighborhood
Row count
Description
least expensive
1763
Houses in a neighborhood with average sale price of under 100,000 USD
less expensive
11746
Houses in a neighborhood with average sale price between 100,000.01 and 250,000 USD
expensive
8301
Houses in a neighborhood with average sale price between 240,000.01 and 500,000 USD
more expensive
1445
Houses in a neighborhood with average sale price between 500,000.01 and 750,000 USD
most expensive
626
Houses in a neighborhood with average sale price of over 750,000.01 USD
Show the code
sum_stat_data %>% dplyr::select(quality_grade) %>%mutate(quality_grade =factor(quality_grade, levels =c("Lowest","Mid","High"))) %>%table() %>%as.data.frame() %>%mutate(Description =c("Houses with a quality grade of \"X\"","Houses with a quality grade of \"X-\", \"A+\", \"A\", or \"A-\"","Houses with a quality grade of \"C\", \"C-\", \"D+\", \"D\", \"D-\", \"E+\", \"E\", or \"E-\"")) %>%kbl(caption ="Table 4: Descriptive statistics of quality grade classes", align =rep("c", 8), col.names =c("Quality Grade", "Row count", "Description")) %>%kable_classic(full_width = F, html_font ="Cambria")
Table 4: Descriptive statistics of quality grade classes
Quality Grade
Row count
Description
Lowest
776
Houses with a quality grade of "X"
Mid
23045
Houses with a quality grade of "X-", "A+", "A", or "A-"
High
60
Houses with a quality grade of "C", "C-", "D+", "D", "D-", "E+", "E", or "E-"
Checking for Collinearity
When using numerous predictors in a model, it is important to inspect for collinearity, or how similar or correlated predictors are to each other. We explored these relationships between features by examining a correlation matrix. As shown in Figure 2, some variables such as number of rooms, number of bedrooms and number of bathrooms have strong collinearity. There is also a high correlation between median household income, spatial lag in price, and percent of population that is white. Taking note of these patterns, we kept our subset of predictors as our next step accounted for collinearity.
It’s also important to contextualize our predictors by looking at the spatial process of these features. It’s very easy to identify the main commercial centers in Philadelphia when mapping out houses’ distance to commercial centers as seen in Figure 3. We see these commercial centers clustered along some main streets in Philadelphia.
Show the code
# define color palettesmypalette1 <-colorRampPalette(c("#bbf0e1","#025e44"))(5)mypalette2 <-colorRampPalette(c("#f0cab9","#6e2e10"))(5)# add geometry back inimp_data_geom <-cbind(imp_data, data_to_dummy %>% dplyr::select(geometry)) %>%st_as_sf() %>%mutate(num_rooms =as.numeric(number_of_rooms))# map of dist_to_commerceggplot() +geom_sf(data = hoods, color ="white", fill ="#f0cab9") +geom_sf(data = imp_data_geom, aes(colour =q5(dist_to_commerce)), show.legend ="point", size = .25) +scale_colour_manual(values = mypalette1,labels=qBr(imp_data_geom,"dist_to_commerce"),name="Distance to Commercial\nCenter (ft)") +labs(title ="Distance to nearest Commercial Center",subtitle="Philadelphia, PA", # make sure to note that NAs were dropped for visualization purposecaption ="Figure 3") +mapTheme()
Show the code
buildtype_data <- data_to_dummy %>%mutate(buildtype =factor(building_type_price_class, levels =c("least expensive","less expensive","expensive","more expensive","most expensive")))# map of build typeggplot() +geom_sf(data = hoods, color ="white", fill ="#f0cab9") +geom_sf(data = buildtype_data, aes(colour = buildtype), show.legend ="point", size = .25) +scale_colour_manual(values = mypalette1,name ="Price Class of\nHouse Type") +labs(title ="Price Class of House Type",subtitle="Philadelphia, PA", # make sure to note that NAs were dropped for visualization purposecaption ="Figure 4") +mapTheme()
Show the code
# map of exterior_conditionggplot() +geom_sf(data = hoods, color ="white", fill ="#f0cab9") +geom_sf(data = imp_data_geom, aes(colour =q5(num_rooms)), show.legend ="point", size = .25) +scale_colour_manual(values = mypalette1,labels=qBr(imp_data_geom,"num_rooms"),name="Number of Rooms") +labs(title ="Number of Rooms",subtitle="Philadelphia, PA", # make sure to note that NAs were dropped for visualization purposecaption ="Figure 5") +mapTheme()
Map of Home Sale Prices in Philadelphia
Show the code
# create price per square foot variableppsqft <- reg_data %>% dplyr::select(sale_price,total_area) %>%mutate(ppsqft_m =ifelse(total_area !=0& total_area !=1,sale_price/total_area, NA)) %>%# there are some places where total_area == 1cbind(data_to_dummy %>%filter(toPredict =="MODELLING") %>% dplyr::select(geometry)) %>%st_as_sf() %>%filter(!is.na(ppsqft_m)) # dropping missing ppsqft_m values for a cleaner mapggplot() +geom_sf(data = hoods, color ="white") +geom_sf(data = ppsqft, aes(colour =q5(ppsqft_m)), show.legend ="point", size = .25) +scale_colour_manual(values = mypalette2,labels=qBr(ppsqft,"ppsqft_m"),name=bquote((USD/ft^2))) +labs(title ="Home Sale Price in Philadelphia",subtitle="Price Per Square Foot", # make sure to note that NAs were dropped for visualization purposecaption ="Figure 6") +mapTheme()
OLS Model
We used an ordinary least squares (OLS) regression model as our predictive model. The components of an OLS model are the dependent variable, independent variable(s), a coefficient for each independent variable, an intercept, and an error value. In our case, home sale price for a the ith house is the dependent variable (\(Y_i\)), or what we’re trying to predict and all of our predictors that we put into the model are independent variables (\(X_{1_i},X_{2_i},...X_{n_i}\); where there are n predictors). The model calculates coefficients (\(\beta_1,\beta_2,...,\beta_n\)) for each independent variable and the intercept (\(\beta_0\)), where the coefficient is defined as the amount of change in sale price we expect for every unit of change in the independent variable that coefficient is paired with. The intercept is the expected sale price given we have no knowledge of any of the predictors. The error value (\(\epsilon_i\)) represents the difference between the actual home sale price and that predicted by the model, and it varies for every prediction.
data.frame(MAE =sprintf('%.2f',round(mae1,3)),MAPE =round(mape1,3)) %>%kbl(caption ="Table 6: Mean Absolute Error (MAE) and Mean Absolute Percentage Error (MAPE) of one test set",align =c("c","c")) %>%kable_classic(full_width = F, html_font ="Cambria") %>%column_spec(1, width ="10em") %>%column_spec(2, width ="10em")
Table 6: Mean Absolute Error (MAE) and Mean Absolute Percentage Error (MAPE) of one test set
MAE
MAPE
67191.43
0.054
Show the code
customSummary <-function(data, lev =NULL, model =NULL) { mape <-mean((abs(data$obs - data$pred)) / data$obs) *100 mae <-mean(abs(data$obs - data$pred)) rmse <-sqrt(mean((data$obs - data$pred)^2)) rsq <-cor(data$obs, data$pred)^2 out <-c(MAE = mae, RMSE = rmse, Rsquared = rsq, MAPE = mape) out}# train_control <- trainControl(method = "cv", # using leave one out instead of regular cv, so gives better estimate of error# number = 10,# summaryFunction = customSummary)train_control <-trainControl(method ="cv", # using leave one out instead of regular cv, so gives better estimate of errornumber =100,summaryFunction = customSummary)model <-train(sale_price ~ .,data = reg_data %>% dplyr::select(-musa_id),trControl = train_control,method ="lm",na.action = na.exclude)print_model <-as.data.frame(print(model))
Linear Regression
23781 samples
31 predictor
No pre-processing
Resampling: Cross-Validated (100 fold)
Summary of sample sizes: 23542, 23543, 23543, 23543, 23544, 23544, ...
Resampling results:
MAE RMSE Rsquared MAPE
68669.9 122346.3 0.7416438 36.7211
Tuning parameter 'intercept' was held constant at a value of TRUE
Show the code
model_resample <- model$resample# summary(model)# coefficient data from summary of model# summ_coeff <- as.data.frame(summary(model)$coefficients)# # model_for_plot <- lm(sale_price ~ ., # data = reg_data)# ols_plot_resid_fit(model_for_plot)data.frame(MeanMAE = print_model$MAE,SDMAE =round(sd(model_resample$MAE),3)) %>%kbl(caption ="Table 7: Results of 100-fold Cross-Validation",align =c("c","c"), col.names =c("Mean MAE", "SD of MAE")) %>%kable_classic(full_width = F, html_font ="Cambria") %>%column_spec(1, width ="10em") %>%column_spec(2, width ="10em")
Table 7: Results of 100-fold Cross-Validation
Mean MAE
SD of MAE
68669.9
6464.582
Show the code
ggplot(data = model$resample) + ggplot2::geom_histogram(aes(x = MAE), fill ="#027a59") +labs(title ="Distribution of MAE",subtitle ="K-Fold Cross Validation; k = 100",x ="Mean Absolute Error",y ="Count",caption ="Figure 7") +plotTheme()
Show the code
phl.test %>% dplyr::select(sale_price.Predict, sale_price) %>%ggplot(aes(sale_price, sale_price.Predict)) + ggplot2::geom_point(color ="#6e2e10", alpha =0.4) +stat_smooth(aes(sale_price, sale_price), method ="lm", se =FALSE, size =1, colour=mypalette1[2]) +stat_smooth(aes(sale_price.Predict, sale_price), method ="lm", se =FALSE, size =1, colour="#027a59") +labs(title="Predicted sale price as a function of observed price",subtitle="Light teal line represents a perfect prediction; Teal line represents prediction",caption ="Figure 8",x ="Observed Sale Price", y ="Predicted Sale Price") +theme_minimal()
Show the code
resids_data <-cbind(reg.training$model, reg.training$residuals, reg.training$fitted.values, phl.training$musa_id) %>%rename(musa_id =`phl.training$musa_id`) %>%left_join(dummied_data %>% dplyr::select(musa_id, geometry), by ="musa_id") %>%mutate(resids =`reg.training$residuals`,pred_value =`reg.training$fitted.values`) %>%st_sf()ggplot() +geom_sf(data = hoods, color ="white", fill ="#f0cab9") +geom_sf(data = resids_data, aes(colour =q5(resids)), show.legend ="point", size = .25) +scale_colour_manual(values = mypalette1,labels=qBr(resids_data,"resids"),name="Residuals") +labs(title ="Map of Residuals for Test Set",subtitle="Philadelphia, PA", # make sure to note that NAs were dropped for visualization purposecaption ="Figure 9") +mapTheme()
Show the code
# adapted from code chunk: autocorrdata_nb <- phl.test %>%left_join(reg_data_geom, by =c("musa_id"="musaID")) %>%st_sf() %>%mutate(nb =st_knn(geometry, k =5),wt =st_weights(nb),.before =1)# calculate global moran's i for sale_priceglbl_moran <-global_moran(data_nb$sale_price, data_nb$nb, data_nb$wt)$`I`# 0.5687531# moran monte carlomoranMC =global_moran_perm(data_nb$sale_price, data_nb$nb, data_nb$wt, alternative ="two.sided", 999)moranMC
Monte-Carlo simulation of Moran I
data: x
weights: listw
number of simulations + 1: 1000
statistic = 0.55464, observed rank = 1000, p-value <
0.00000000000000022
alternative hypothesis: two.sided
Show the code
moranMCres = moranMC$res |>as.data.frame() %>%rename(moranMCres =1)ggplot(moranMCres) +geom_histogram(aes(moranMCres), bins =100, fill ="#027a59") +geom_vline(xintercept = glbl_moran, col ="#6e2e10") +labs(title ="Observed and Permuted Moran's I",subtitle ="Brown line indicates the observed Moran's I",x ="Moran's I",y ="Count",caption ="Figure 10") +plotTheme()
Show the code
# this can be made prettierspdep::moran.plot(data_nb$sale_price, nb2listw(data_nb$nb),xlab ="Sale Price", ylab ="Spatial Lag")
Show the code
# add caption that says Figure 11
TO DO: Provide a map of your predicted values for where ‘toPredict’ is both “MODELLING” and “CHALLENGE”
Show the code
inTrain2 <-createDataPartition(y =paste(imp_data$fireplaces, imp_data$number_of_rooms, imp_data$exterior_condition), p = .60, list =FALSE) # p = .60 = split to 60/40 ratio for train/testphl.training2 <- imp_data[inTrain2,] # defining training setphl.test2 <- imp_data[-inTrain2,] # defining testing setreg.training2 <-lm(sale_price ~ ., data =as.data.frame(phl.training2) %>% dplyr::select(-c(musa_id, to_predict)))# plot(reg.training2) # useful plots to look at model fitphl.test2 <- phl.test2 %>%mutate(sale_price.Predict =predict(reg.training2, phl.test2), # ^^ actually adding the predicted values from model defined in reg.training2sale_price.Error = sale_price.Predict - sale_price,sale_price.AbsError =abs(sale_price.Predict - sale_price),sale_price.APE = (abs(sale_price.Predict - sale_price)) / sale_price.Predict) %>%left_join(dummied_data %>% dplyr::select(musa_id, geometry), by ="musa_id") %>%st_sf()mape2 <-mean(phl.test2$sale_price.APE) # 0.3139695mae2 <-mean(phl.test2$sale_price.AbsError) # 69252.47summ_reg.training2 <-summary(reg.training2)summ_coeff2 <-as.data.frame(summary(reg.training2)$coefficients)ggplot() +geom_sf(data = hoods, color ="white", fill ="#f0cab9") +geom_sf(data = phl.test2, aes(colour =q5(sale_price.Predict)), show.legend ="point", size = .25) +scale_colour_manual(values = mypalette1,labels=qBr(phl.test2,"sale_price.Predict"),name="Predicted Sale Price") +labs(title ="Map of Predicted Sale Price for Test Set",subtitle="Test Set includes both modelling and challenge data; Philadelphia, PA", caption ="Figure 12") +mapTheme()
In comparison to the map of observed sale prices (figure 6), we can see in this map of the predicted sale price (including both modelling and challenge data; figure 12) that our model generally predicts prices well, but overpredicts sale price in the northeast region as well as the northwestern tip of Philadelphia.
TO DO: Using the test set predictions, provide a map of mean absolute percentage error (MAPE) by neighborhood. (code chunk reids is close to what we want, but use MAPE instead of resids)
Show the code
# code adapted from resids code chunk + book# by "test set" do they mean the new test set with both modelling and challenge data?# proceeding with old test set (only modelling) for nowmape_data <- phl.test2 %>% dplyr::select(sale_price,sale_price.Predict,sale_price.APE,musa_id) %>%left_join(dummied_data %>% dplyr::select(musa_id, geometry), by ="musa_id") %>%st_sf()mape_data <-st_join(mape_data, hoods)mape_x_hood <-st_drop_geometry(mape_data) %>%group_by(mapname) %>%summarize(mean.MAPE =mean(sale_price.APE, na.rm = T),mean.price =mean(sale_price, na.rm = T),mean.pred.price =mean(sale_price.Predict, na.rm = T)) mape_x_hood %>%ungroup() %>%left_join(hoods) %>%st_sf() %>%ggplot() +geom_sf(aes(fill = mean.MAPE)) +# geom_sf(data = bothRegressions, colour = "black", size = .5) +scale_fill_gradient(low = mypalette2[1], high = mypalette2[5],name ="MAPE") +labs(title ="Mean test set MAPE by neighborhood",caption ="Figure 13") +mapTheme()
Show the code
# mape_x_hood <- mape_data %>%# group_by(mapname) %>%# summarize(avg_mape = mean(sale_price.APE),# avg_price = mean()) %>%# arrange(desc(avg_mape)) %>%# select(avg_mape, mapname) %>%# st_drop_geometry() %>%# left_join(., hoods, by = "mapname") %>%# st_sf()# # # redefine qBr to show decimal points# qBr_3 <- function(df, variable, rnd) {# if (missing(rnd)) {# as.character(quantile(round(df[[variable]],3),# c(.01,.2,.4,.6,.8), na.rm=T))# } else if (rnd == FALSE | rnd == F) {# as.character(formatC(quantile(df[[variable]],# c(.01,.2,.4,.6,.8), na.rm=T),# digits = 3))# }# }# # ggplot() +# # geom_sf(data = hoods, color = "white", fill = "#f0cab9") + # geom_sf(data = mape_x_hood, aes(fill = q5(avg_mape)), # show.legend = "point", size = .25) +# scale_fill_manual(values = mypalette1,# labels=qBr_3(mape_x_hood,"avg_mape"),# name="MAPE") +# labs(title = "Map of MAPE by Neighborhood",# subtitle="Test Set includes both modelling and challenge data; Philadelphia, PA", # caption = "Figure 13") +# mapTheme()
TO DO: Provide a scatterplot plot of MAPE by neighborhood as a function of mean price by neighborhood. aki: are they talking about predicted price or actual price? actual price plotted for now
Show the code
st_drop_geometry(mape_x_hood) %>%ggplot(aes(mean.price, mean.MAPE)) + ggplot2::geom_point(size = .5, color ="#6e2e10", alpha =0.4) +geom_smooth(method ="lm", se=F, colour ="#027a59") +labs(title ="Mean MAPE by neighborhood as a function of mean price by neighborhood",y ="Mean MAPE by neighborhood", x ="Mean price by neighborhood",caption ="Figure 14") +plotTheme()
Show the code
# possibly add neighborhood names of extreme outliers
TO DO: Using tidycensus, split your study area into two groups (perhaps by race or income) and test your model’s generalizability. Is your model generalizable?
st_join(mape_data, split_data) %>%filter(!is.na(incomeContext)) %>%group_by(incomeContext) %>%summarize(mean.MAPE = scales::percent(mean(sale_price.APE, na.rm = T))) %>%st_drop_geometry() %>%# spread(incomeContext, mean.MAPE) %>%kbl(caption ="Table 8: Test set MAPE by neighborhood income context", align =rep("c", 8),col.names =c("Income Context", "Average MAPE")) %>%# footnote("70 rows of data were missing median household income and are therefore excluded from this table.") %>% kable_classic(full_width = F, html_font ="Cambria") %>%column_spec(1, width ="12em") %>%column_spec(2, width ="12em")
Table 8: Test set MAPE by neighborhood income context
Income Context
Average MAPE
High Income
23%
Low Income
56%
Show the code
st_join(mape_data, split_data) %>%group_by(raceContext) %>%summarize(mean.MAPE = scales::percent(mean(sale_price.APE, na.rm = T))) %>%st_drop_geometry() %>%kbl(caption ="Table 9: Test set MAPE by neighborhood race context", align =rep("c", 8),col.names =c("Race Context", "Average MAPE")) %>%# footnote("70 rows of data were missing median household income and are therefore excluded from this table.") %>% kable_classic(full_width = F, html_font ="Cambria") %>%column_spec(1, width ="12em") %>%column_spec(2, width ="12em")
Table 9: Test set MAPE by neighborhood race context
Race Context
Average MAPE
Majority Non-White
58%
Majority White
22%
Show the code
# data_nb <- data %>% # mutate(nb = st_knn(geometry, k = 5),# wt = st_weights(nb),# .before = 1)# # # # calculate global moran's i for sale_price# global_moran(data_nb$sale_price, data_nb$nb, data_nb$wt)$`I`# # # moran monte carlo# moranMC = global_moran_perm(data_nb$sale_price, data_nb$nb, data_nb$wt, alternative = "two.sided", 999)# moranMC# # moranMCres = moranMC$res |># as.data.frame()# # colnames(moranMCres) = "col1"# # ggplot(moranMCres) +# geom_histogram(aes(col1), bins = 100) +# geom_vline(xintercept = moranMC$statistic, col = "red") +# labs(title = "Histogram of MoranMCres",# x = "moranMCres",# y = "Frequency",# caption = "Figure xx") +# theme_minimal() +# theme(plot.title = element_text(hjust = 0.5))# # # spdep::moran.plot(data_nb$sale_price, nb2listw(data_nb$nb),# xlab = "Sale Price", # ylab = "Spatial Lag")
Is this an effective model? What were some of the more interesting variables? How much of the variation in prices could you predict? Describe the more important features? Describe the error in your predictions? According to your maps, could you account the spatial variation in prices? Where did the model predict particularly well? Poorly? Why do you think this might be?
Accuracy and Generalizability
Conclusion
Would you recommend your model to Zillow? Why or why not? How might you improve this model?
Footnotes
Michael Neal, Sarah Strochak, Linna Zhu, and Caitlin Young, “How Automated valuation models can disproportionately affect majority-Black neighborhoods,” Urban Institute, December, 2020.↩︎
Alan S. Dornfest, Chair Bill Marchand, Doug Warr, August Dettbarn, Wayne Forde, Joshua Myers, and Carol Neihardt, Standard on automated valuation models (AVMs), International Association of Assessing Officers, July, 2018.↩︎
